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Abstract 

We begin an investigation of inhomogeneous structures in holographic su- 
perfluids. As a first example, we study domain wall like defects in the 3+1 
dimensional Einstein-Maxwell-Higgs theory, which was developed as a dual 
model for a holographic superconductor. In [1], we reported on such "dark 
solitons" in holographic superfluids. In this work, we present an extensive nu- 
merical study of their properties, working in the probe limit. We construct 
dark solitons for two possible condensing operators, and find that both of them 
share common features with their standard superfluid counterparts. However, 
both are characterized by two distinct coherence length scales (one for order 
parameter, one for charge condensate). We study the relative charge depletion 
factor and find that solitons in the two different condensates have very distinct 
depletion characteristics. We also study quasiparticle excitations above the 
holographic superfluid, and find that the scale of the excitations is comparable 
to the soliton coherence length scales. 



* ville . keranen@helsinki . fi 

' esko . keski- vakkuri® helsinki . fi 

* sean. nowling@helsinki . fi 

§ yogendran. kalpat @helsinki . fi 

1 



1 Introduction 



The AdS/CFT correspondence [2-4] - or gauge-gravity duality, as its more phe- 
nomenological incarnation is called - is being used as a tool to study an increasing 
variety of strongly interacting systems. When the system is also " close" to being con- 
formal, it is natural to try to construct a holographic gravity dual model. The main 
practical advantage of holography is that it allows a new, computationally simpler 
framework to explore quantum systems which are outside the reach of conventional 
perturbative techniques. The situation of interest in this document will be a system 
in which a U(l) symmetry is spontaneously broken, as is relevant for superfluidity 
and superconductivity. 

Superfluidity may occur in a system at low temperatures when the ground state 
becomes occupied by a macroscopic number of particles. This is known to occur in 
both interacting bosonic systems (Bose-Einstein condensates, BEC) and interacting 
fermions with Bardeen-Cooper-Schrieffer pairing (BCS). In fact, cold atomic systems 
can display both BEC-like and BCS-like superfluidity and even a smooth transition 
between them. 

It is not at all obvious what kinds of superfluids may be realized holographically. A 
prototype example is the recently constructed gravitational dual theory [6], building 
on the model [7], for the purpose or modeling BCS superconductors [6,8] or relativis- 
ts superfluids with a spontaneously broken global U(l) symmetry [9]. The model 
was also extended to an effective theory for some properties of a class of quantum 
Hall fluids [10]. Typically, in systems that are dual to AdS models, one finds both 
fermionic and bosonic excitations. Thus, one might hope to be able to holographi- 
cally model both BEC-like and BCS-like superfluidity, as well as more complicated 
fluids with mixed behavior. For closer contact with real world cold atom systems, a 
required ingredient would also be manifest non-relativistic symmetry, for which the 
first holographic models were studied in [11, 12]. In this paper, and in the associated 
sequence of works to follow, we will focus on studying extended configurations in the 
model [6], and we will find a richer variety of possibilities than previously expected. 

Hydrodynamic properties of the relativistic holographic superfluid have been ex- 
plored in [13,14] - the latter in particular argues that the superfluid does not obey the 
Landau relationship between sound velocities. In [15], masses of quasiparticles were 
determined analytically near T c and substantially extended in [16] (to other values of 
the mass of the bulk scalar field as well). 

Another basic fact is that conventional superfluids are known to support long lived 
spatially inhomogeneous configurations with nontrivial topology. An experimentally 
and theoretically interesting class of defects, found in superfluids, are domain wall 
like defects called dark solitons. A dark soliton is an interface of reduced (charge) 
density between two superfluid phases, with the order parameter changing sign across 
the interface. Holographic counterparts were recently reported by us in [1]. Apart 
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from curiosity, there are several reasons to study such solitons, a primary one being 
to explore their properties at strong coupling, a difficult problem using traditional 
methods. Although the gravity dual has different symmetries and microscopic prop- 
erties from the real life cold atomic systems, there may be some universal features. At 
the very least it provides a "spherical cow" example allowing one to study properties 
of a class of dark solitons at strong coupling. 

The holographic dark solitons are also interesting solutions in their own right - 
they turn out to reveal important information about the different ways of condensing 
operators in [6] and about the nature of the holographic superfluid. First of all, we 
find that the solitons have two distinct length scales: one for the charge density and 
one for the order parameter. This behavior is in contrast to a single scale found in 
simple solutions to the Gross-Pitaevskii equation for superfluids. Furthermore, by 
studying the degree of charge density depletion we discover that the two different 
choices for the condensing operator, discussed in [6], lead to different properties for 
the associated superfluids. We will also compare the length scales as set by these 
solitonic configurations with the microscopic masses of quasiparticles. 

This document is organized as follows. We begin with a brief discussion of the 
holographic description of a superfluid. The configurations of interest are obtained 
by solving a system of partial differential equations in AdS space. These equations 
and the boundary conditions - which sustain the soliton - are presented in Section [2j 
However, the equations seem intractable analytically. We therefore have to numeri- 
cally solve the equations. The discretization of the equations and the method we use 
to solve the equations are presented in Section [31 In order to identify properties of 
interest to be obtained from the numerical solutions, we will briefly discuss the Gross- 
Pitaevskii equations in Section HI We will then present our results for the holographic 
dark solitons and conclude with a discussion of the results and future directions. 



2 Holographic Description 

An explicit holographic modeling of a superfluid system was constructed in [6] fol- 
lowing the ideas in [7] (closely related ideas also appeared in [17]). One considers 
a system consisting of a complex scalar field \l/ interacting with a U(l) gauge field 
Ffj, u — d^A u — d u A^ in 3+1 dimensions in the presence of gravity with a cosmological 
constant L. The action for this system is 



^{R - 12A) + \ {--.F^ - m 2 M - D^D^ 
q \ 4 

where the covariant derivative = (<9 M — iA^)^ and q may be identified with 

the charge of the scalar fielcEJ. Focusing on just these degrees of freedom (rather 



1 We have rescaled the fields such that q only appears in front of the matter action. 
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than keeping all the fields required for a UV complete gravity theory), amounts to 
concentrating on the most relevant operators of the field theory. One could include 
other fields (say, a fermion field which would be dual to fermionic quasiparticles). 
In keeping with the usual terminology, we will refer to the fields of the gravitational 
theory as the "bulk" fields in what follows. 

The AdS / CFT correspondence [2] requires us to find solutions to the equations of 
motion of the 4D-gravitating system, such that the metric is asymptotically AdS. The 
properties of the dual field theory are then "read off" from the asymptotic behavior 
of the various fields of this gravitating system ( [3,4], see the review [5]). 

In our case, since we are trying to study a spatially inhomogeneous configuration, 
solving the full gravitating system is a difficult problem. Therefore, we consider a 
limit in which the energy density in the complex scalar field and the Maxwell field 
(in 4D) is scaled to zero. This may be achieved by taking — — > keeping the field 
values finite. In this limit (termed the probe limit in the AdS literature) the 

metric is determined, independent of the matter fields, to be a 4-D planar AdS black 
hole with a metric 

, 9 r9 . fdt 2 dz 2 dx 2 . „. . _ ( z \ 3 

& = l\->— + - 2 + -), /w = i-(-J • (2) 

The Hawking temperature of the black hole, T H = is identified as the equilibrium 
temperature of the dual field theory. 

The complex scalar field and gauge field then propagate on this background and 
their dynamics is determined by an action 

S = J d^xdz^g (^—F^F^ - m 2 ^ - D^D^ (3) 

where the covariant derivative = (<9 M — iA^ and the various indices are con- 
tracted using the metric given above. 
If we rescale 

A * 

(z, x) -> z T (z, x) A -> — V -> — , (4) 

Zj* Zt 

then the new co-ordinates z,x are dimensionless as are the new fields A. Further, 
all dependence on zt is also removed. 

The equations of motion for this system are 

= -^^(v^OTj-m^ (5) 

= ^7=^ {V^gF^) + i(yd v y - $>d v y) + 2A V ^. (6) 

In this work we will primarily be interested in studying static kink solutions in the field 
theory and therefore the fields only have dependence on z and one spatial variable, 
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x. In the gauge where A z = 0, the A z equation is satisfied if \I/ is taken to be real 
and A x = 0. We will assume translation invariance in the ^-direction and set A y = 0. 
Then the only nonzero component of the gauge field is A = A(z,x) 

It is important to note that, after gauge fixing, there is a residual Z2 symmetry 
under which the scalar field changes sign. When studying dark soliton solutions to the 
equations of motion, the scalar field changes sign as a function of x, lim^-too = 
±|\]/(oo)|. One might think that this sign flip may also be gauged away, but this would 
require a gauge transformation that is singular as the scalar field crosses through zero, 
and hence is not an allowed gauge transformation. 

It is convenient to redefine the field slightly, \I/ = zR/y/2, and rewrite equations 
of motion 

fA" + d 2 x A -R 2 A = (8) 
2.1 Bulk- boundary dictionary 

Typically, in AdS spacetimes, the solutions to the equations of motion may be segre- 
gated into " normalizable" and "non-normalizable" parts determined by their leading 
z behavior as one approaches the boundary. According to the AdS/CFT dictio- 
nary [3-5] the two have different interpretations in the dual field theory. The bound- 
ary values of the non-normalizable modes are interpreted as sources for the operators 
of the dual field theory. The normalizable modes are then identified with the vacuum 
expectation values of the corresponding operators sourced by the non-normalizable 
modes [18]. The scaling behavior of the dual operator is also fixed and may be ex- 
tracted from the asymptotic behavior of the normalizable modes as functions of z. 
This is explicitly realized by identifying the generating functional of connected cor- 
relation functions, W[J], of the strongly coupled field theory with the on-shell action 
for the bulk gravity theory (the bulk fields have boundary values (pg = J. ( [3,4]) 

S Gl&y .,o.s.\ <t>d=J = W QFT [J]- (9) 
One simple way to see the relationship between normalizable modes and expectation 
values is to consider the change in the action of the gravity theory caused by infinites- 
imal changes in the "non-normalizable" gravity modes (field theory sources). Upon 
partial integration and the inclusion of a suitable counterterm, the variation of the 
action ([I]) in terms of the fields A, becomes 

5S = - J d A xd^y/gd^) + ^(y/gdA F»°) (10) 

+ J d 4 x5V + Ua v £{A) 

+ terms higher order in the variations. 
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We have grouped the equations of motions in ((Sj) as S(^f) and S(A) for brevity. 
Rewriting the surface terms as in the AdS prescription gives us 



SS = - J d 3 x[(5^^d z ^) + (y/g5A F z0 )} 
+ J d 4 x8V 5(f) + -8 A £{A) 
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+ terms higher order in the variations. 

contributions only from the first line in the above. 

As in [6], we shall take the scalar field to satisfy m 2 = —2/L 2 . For this mass value, 
there are two quantizations for scalar fields in any asymptotically AdS4 spacetime [19]. 
The two quantizations correspond to exchanging the role of source and expectation 
value in the dual field theory. 

More explicitly, it can be easily seen that, close to the boundary at z = 0, solutions 
to the equations of motion must behave as 

R ~ R w + zR® + A ~ A^ + zA^ + ... (12) 

in an expansion of the Frobenius type along the z-direction. We can reg ard R^ as 
the source for a charged, dimension 2 operator 

rM = ztJ® R {2) = z 2 (0 2 ) Af = z T p A^ = z 2 p. (13) 

where (O2) is the charged operator, /1 is the chemical potential, and p is the charge 
density (all in the dual field theory). One could motivate the identification of the 
asymptotic value A^ = z T p by noting that, in Euclidean space, a vev of A is 
naturally interpreted as a chemical potential (since it minimally couples to a conserved 
charge) . 

In the second quantization scheme, we identify w 2 ^ 1 as the source for a charged, 
dimension 1 operator 

RV> = *r(O x ) RW=4JV> Af=z T p Af = z\p. (14) 

(we will use ip to denote the condensate in what follows). The various prefactors 
come from the rescaling of z by z? so that the fields in the bulk are dimensionless. 
Solutions to the gravitational equations of motion with nonzero non-normalizable 
components are interpreted in the dual field theory as deforming the Hamiltonian, 
57i ~ J dPxJ^Ow and would hence change the ground state of the theory. 

To obtain a unique solution to the bulk equations of motion, we need to impose 
boundary and regularity conditions. As in [6], we want to study the spontaneous 
breaking of the U(l) symmetry unperturbed by any sources, so we will impose bound- 
ary conditions such that in each quantization the source terms are turned off. 
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In addition, we impose the following regularity conditions at the horizon (we 
assume that R is regular at the horizon) 

f'R' - zR + d 2 x R = A (z = l,x) = (15) 

(the condition on the gauge field has been argued to be a regularity condition [7]). Be- 
cause the differential equations are elliptic ([7]) outside the horizon, Dirichlet /Neumann 
boundary conditions along the boundary (at z = 0) and regularity (at the horizon) 
are sufficient to fix a unique solution (a useful analogy is with Poisson equations). 



2.2 Review of earlier results 

One of the features of the AdS/CFT correspondence is that global symmetries of the 
field theory system appear gauged in the gravitational theory. Given this, the U(l) 
gauge symmetry in the gravitational Lagrangian is interpreted as corresponding 
to a global U(l) symmetry of the dual field theory system. And hence, the bulk 
charged scalar field \I/ maps into a field theory operator Oi that is "charged" under 
the global U(l). We can now study this system to see if it has ground states in 
which these charged operators pick up nonzero vacuum expectation values (vevs) in 
the limit where all external sources for charged operators are turned off. 

In the works of [6], it was shown that for fi > fi c (or equivalently small T, since 
there is only one independent parameter p, = J^), one can find nontrivial solutions for 
the scalar field equations displaying spontaneous symmetry breaking in the dual field 
theory. Thus, for low enough temperature (or large enough /i) we obtain superfluidity 
(strictly, we have only argued for a condensate, that this is indeed a superfluid has also 
been established by showing the existence of a hydrodynamic mode [9]). The critical 
value of [L (or T) depends on whether we consider Rp^ or R^ as our normalizable 
mode (or equivalently, as the vev of the order parameter). 

The graph of typical solutions look as in Fig. [1] wherein we plot the two solutions 
obtained by condensing either operator 0\ (dotted line) or O2 (dashed line) We also 
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Figure 1: Scalar field and Gauge field profile: The dotted line refers to 0\ and the 
dashed line to 2 type condensate 



note that at any temperature (or chemical potential) we have a solution where the 
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condensate is absent A(z) = p,(l — z) and ^ = 0. Such a solution is energetically 
disfavored at low temperatures (T < T c or /x > /i c ). [9] 

The features in the bulk gauge field are highlighted better if we graph a subtracted 
potential AA(z) = A(z)/A(0) — (1 — z) which may be interpreted as charge density 
over and above the zero condensate value. It is seen in Fig. [1] that the bulk gauge 
field profile is quite similar in the two cases. One feature of potential interest is that 
the presence of the condensate changes the charge density by little (vide the range of 
the vertical axis), at least for small chemical potentials. 

We also plot the charge density and condensate as a function of the holographic 
direction z. One of the interesting differences between the two types of solutions is 
that when O2 condenses, the condensate always vanishes at a point z = z c . Beyond 
this z c value the condensates is always monotonically decreasing as one approaches 
the horizon. For 0\ condensates, the condensate is monotonic as a function of z. 
The Landau-Ginzburg approach suggests that the square of the condensate should be 
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Figure 2: The condensate and charge density profile: The dotted line refers to 0\ 
and the dashed line to 2 type condensate 



function of (/i — /i c ) ( |26|) . As can be seen from Fig. [31 for small values of , the 
curve may be approximated by a straight line. However, for somewhat larger values, 
the slope changes - and the fitting functions require higher powers. 

The equation of state for this system - namely the graph of the charge density 
as a function of the chemical potential for fixed temperature looks as in the Fig. |H 
Close to the phase transition, the charge density and chemical potential are linearly 
related, but near T = (large /i), this is no longer true. 

One could expect that the condensate should be a linear function of the charge 
density - if one imagines that the microscopic system consists of fermions, then both 
the charge density and the condensate are proportional to the vev of ipip - where ip is 
the fermion many body wavefunction. The condensate is indeed a linear function of 
the charge density at large values of the chemical potential (i.e., close to T = 0) but 
near the phase transition it seems that increasing the condensate does not increase 
the charge density by a considerable amount. 
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Figure 3: The condensate as a function of the chemical potential: Upper curve is O2, 
lower 0\. 
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Figure 4: The charge density vs. chemical potential: Upper curve is O2, lower 0\. 

3 Numerical Methods for Dark solitons 

The bulk equations of motion (J7J) seem to be intractable even in the homogeneous case. 
Therefore, we resort to numerical methods to solve them. Since we are interested in 
inhomogeneous solutions, shooting methods become unwieldy (in this regard, see [33] 
for a tour de force of numerical work in a closely related context). Therefore, because 
the differential equations are elliptic outside the horizon, it is natural to use relaxation 
methods to find solutions. 

More precisely, to simulate the differential equations, we place the system in a 
large box of size 1 x L x in the z and x directions respectively. We discretize this box 
using a lattice, and then solve a discretized version of the differential equation on 
the lattice using a Gauss-Seidel routine. The Gauss-Seidel strategy is to start with 
a seed configuration which obeys appropriate boundary conditions. We can then use 
a discrete representation of the differential equations to "relax" this configuration 
towards a solution. If the lattice is sufficiently fine, the elliptic nature of the problem 
implies that the seed configuration will flow to an exact solution to the differential 
equation eventually. 
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Figure 5: The condensate as a function of the charge density (upper curve is 0\ 
condensate) 



We set our conventions by labeling the z and x positions with i e [0, M] and 
j G [0, N] respectively and use step sizes h z and h x . We will denote the lattice fields 
as (R, A) -> (Ri,j, Aij). 

The discretization of these equations (|7j) has to be done separately for the interior 
of the lattice and the edges to at least second order in the lattice spacings. In the 
interior, when i ^ 1, M j ^ 1,N, we can use the center difference formulae for the 
derivatives, 



df(x) 



dx 
d 2 f(x) 



2h a 



0{h\ 



dx 1 



f(x i+1 )-2f(x i ) + f(x i - 1 ) 
hi 



+ 0(h 2 x ) 



(16) 



and therefore we get the following algebraic equations 



ft 



(( 



h 2 , 2h z 



)R+0 + 



■J_ 

-hi 



JL 

2h z 



)Rn + 



Rq++Rq- 
hl 



00 



.4, 



^-^ + 4 + 1) 

- ) + ±(A 0+ + A _ 



(17) 



'.hi 



( /?2 il + I 



The above equations are to be understood as follows. At a site the value of 
the scalar field Rij which we label as i?oo is determined in terms of its neighbors 
Ri±i,j — R±o etc by the above equations. In addition, z refers to the lattice value Zi, 
and / refers to f(zi) and /' must be replaced by the center difference formula at the 
lattice site. 

At the spatial edges (in the bulk) when x = ±L x /2 and z 7^ 0,1, we simply 
use one-sided representations of the finite difference derivatives and again impose the 
equations of motion. We have also checked that one could also employ Neumann 
boundary conditions at x = without changing the numerically determined solu- 
tions, as long as L is sufficiently large. 
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At the horizon, we impose regularity conditions (|T5|) appropriately discretized. At 
the boundary of AdS space, we have to numerically impose the boundary conditions 



(D\ case : A(z — 0, x) = \i — constant, ~g~( z = 0, x) = (19) 
2 case : A(z — 0, x) = \i — constant, R{z — 0, x) = 0. 

By cycling through the lattice, imposing either the equations of motion, regularity, 
or boundary conditions depending upon the position, an initial seed configuration 
relaxes to a solution to the algebraic equations satisfying the appropriate boundary 
conditions. 

As a first check of our algorithm, we checked that we can obtain the symmetry 
breaking solutions obtained by [6] that correspond to having a homogeneous phase 
without any dark solitons. For the O2 condensate for instance, we set the chemical 
potential at the boundary to a fixed value, and -R(0, x) = 0, and allow fairly arbi- 
trary seed field values in the bulk of the lattice (as well as at the horizon). Upon 
subsequent iteration, we obtain a translationally invariant solutions that match the 
numerical solutions obtained by solving the corresponding one dimensional problem 
using Mathematica's NDSolve to remarkable accuracy even on fairly modest lattices. 
For the 0\ condensate we enforced d z R(0, x) = instead. Again, after iterating one 
finds spatially homogeneous solutions which agree with [6]. 

3.1 Constructing Dark Solitons 

The dark solitons were constructed, numerically, as follows. We first chose a seed 
configuration that asymptotes to solutions in [6] far away from the interface. The 
initial configuration is assumed to be odd in the x-direction. We then numerically 
iterate the seed configurations until it relaxes to a stable configuration (we do this on 
several lattice sizes for each soliton). In actual calculations, we are forced to use a 
cutoff at a distance e = 10~ 10 from the horizon at z = 1. We have checked that this 
value of e does not affect the results quoted. 

A typical solution so obtained is shown in Fig. E] for the 0\ case. In the gauge 
field plot, we have subtracted a linear part — z). Since this is an 0\ condensate, 
the boundary values of the scalar field are non-vanishing. A corresponding typical 
solution for the O2 case is shown in Fig. [7J In this case the boundary value of the 
scalar field is zero and the derivative of the scalar field, dzR, is identified with the 
condensate, magnitude of 

Further, starting with one such solution, we can perturb the solution by an arbi- 
trary deformation that preserves the boundary values. We have seen that the per- 
turbed configuration relaxes back to the original starting solution rapidly (in iteration 
time). This is an indirect argument for the stability of the solution as well. 
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Figure 6: 0\ case: The bulk solution (Scalar field on the left and subtracted gauge 
field on the right). 




Figure 7: O2 case: The bulk solution (Scalar field on the left and subtracted gauge 
field on the right). 

3.2 Numerical Errors 

There are two main sources of errors in the relaxation method described above. First, 
there is the discretization error caused by using a lattice description. This may be 
countered in two obvious ways. We could use discrete representations of derivatives 
that are of higher order in h x and h z and/or use a finer lattice. We however resort 
to a version of the multigrid method. That is to say, we first solve the numerical 
system on a coarse lattice and then use this approximate solution to seed the starting 
configuration on a finer lattice. This procedure converges much faster than naively 
solving on the larger lattice. The second source of error is the degree to which we 
solve the algebraic equations obtained after discretizing the differential equations. 
The most obvious way to minimize the algebraic error is to simply let the system 
relax over more iterations. 

As a measure of the error in our numerical solutions, we define a quantity which 
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we call equation error, 8. This is the extent to which the numerically determined 
solutions do not solve the equations of motion schematically written as 

E.O.M(R num ,A num )=S. (20) 

To evaluate £, we can again discretize the equations of motion by now using n-th order 
representations for the various derivatives and evaluate S num (n) = E.O.M n (R num , A num ). 
For a given numerical solution, we find that S num (n) hardly changes for n > 5. We 
may therefore use this numerically determined equation error (at large derivative 
ordern > 5 ) as an approximation of the exact error in the numerical solution 

• (21) 

large n 

In order to reduce the discretization error we may pass numerically obtained solutions 
obtained on coarse lattices on to finer lattices and repeat the process. In this manner 
it is possible to systematically reduce the equation error, £. In principle this process 
could be continued indefinitely, although in reality we are limited to lattices of size 
less than 513 x 513 sites when using Mathematica on a desktop computer. 



Errors 

Eqn Error 




Figure 8: £dis(7) as a function of lattice size. The dashed line is the maximum value 
of the R equation of motion evaluated on our numerical solution and the solid line is 
the corresponding value for the gauge field 



3.2.1 Uncertainty Analysis 

Although we have introduced a measure of numerical error in the solutions, it would 
be more insightful to have a measure of the error in physical quantities. One way 
to do this comes from (fTTl) which establishes the relationship between normalizable 
modes and operator expectation values. Specifically, the error in the condensate and 
charge density may be estimated by noting that since the equations of motion do not 
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actually evaluate to zero, 



A«0 2 » 
A«P» 



d x 



5R 



6RV* 



r d 4 x 


5R 


1 f(z) 


5AM 



.(£)! + 



.(£)! + 



5A 



6RW 
5A 



\£num(.A^) | 



SAW 



(22) 
(23) 



for the O2 condensate. Similar expressions may be obtained for the 0\ condensate. 

As a final input we need to evaluate the functional derivative of bulk fields with 
respect to their boundary values. We will approximate the functional derivatives by 
their values for the homogeneous solutions (as determined numerically using Mathe- 
matica's NDSolve routine). That is to say, we determine the change 5R(z) when we 
change the boundary values by SR^\ We have checked that this latter approxima- 
tion is a good estimate of the actual functional derivative jj^j even if we use the full 
inhomogeneous equations. 

Such as they are, these error estimates do not represent upper bounds on the 
total error on spatially inhomogeneous solutions, but we do believe that they are a 
descriptive of the error in our numerical solutions. 



4 Dark Solitons 

In order to identify quantities of interest to be extracted from our numerical solutions, 
we will attempt to compare our holographic solitons with the soliton solutions of 
the Gross-Pitaevskii(GP) equation (the latter gives a coarse-grained description of a 
superfluid valid at long wavelengths). However, we emphasize that we will only be 
using the GP equation as a guide - the results indicate that the holographic solitons 
are not the same as the solution of the GP equation. We also note that the description 
of a superfluid system in terms of the GP equation alone is inadequate to capture 
the effects of significant density depletion in the superfluid (perhaps due to strong 
interaction effects) and also the effects of having a significant non-condensate fraction 
(at say high temperatures). In our results, we find that both these contributions seem 
to be present and significant in the holographic superfluids. 

The order parameter of superfluidity is a complex scalar field (also called the 
condensate wave function)- for which one can write down equation of motion called 
the Gross-Pitaevskii equation 

- ^~d 2 i> + (V- + g^\ 2 = (24) 
ZrriB 

(we have dropped the time dependence in the above since we are only interested in 
static phenomena in this work). 
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The dark soliton, in this language is a spatially varying solution of the GP equation 
which interpolates between the potential minima 



ip{x 



oo 



A 



ij)(x — > — oo) 



-A 



(25) 



If we further assume translational invariance in the y-direction, the GP equations 
(1241) has a well known exact solution 



V 



9 



tanh(x/£), 



(26) 



which interpolates between the two minima at ip — ± A = y . The correlation 
length (or the healing length) £, can then be written in a useful form in terms of the 
parameters of the GP equation as 



e 



i 



i 



2gA 2 (V-fi)m B ' 



(27) 



Although the GP equation is really relevant for a nonrelativistic system, it will be 
interesting to test these dependences of the coherence on the magnitude of the con- 
densate and on the chemical potential. Finally, because the charge density is simply 
related to the order parameter p ~ \ip\ 2 , the solution (|26|) has vanishing charge density 
at its core. 

4.1 Holographic results 




Figure 9: 0\ Condensate and charge density as a function of x 



From the numerical solutions fl6|7|) obtained by solving the bulk equations - we 
can extract the boundary profiles of the charge density and the condensate R^ 1 ' 2 ) 
respectively. The plots as well as expectations from the GP equation suggest that the 
condensate can be fitted by a tanh(|) profile and hence plausibly, the charge density 
can be fitted by a sech 2 (^-) profile. These data points along with the best fit curve are 
shown in Fig. [9] for a typical 0\ type of condensate. The results of a similar analysis 
for the 2 system is shown in Fig. [TUJ The profiles of the boundary observables are 
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Figure 10: O2 Condensate and charge density as a function of x. 



equally well fit for both 0\ and O2 cases. Using a least square fit, we then extract 
the coherence length £ from the condensate profile and £ g from the charge density. It 
must be noted here that the fit value of £ and £ 9 are quite sensitive to the accuracy of 
the solutions. We have ensured that the equations are solved to an accuracy of 10~ 4 
- that is to say, the maximum absolute value of the right hand side of the differential 
equation is less that 1CT 4 . 

Following our earlier discussion of the GP-equation ff271) . we plot the behaviour 
of the coherence length as a function of the condensate in Fig. [TTJ One of the first 
surprises is that £ is a linear function of the inverse condensate over a large range of 
chemical potential, although there is no obvious reason for it in this context. 




1 2 3 4 5 f t 



Figure 11: O2 Coherence length versus inverse condensate. 



4.2 Two Length Scales 

Secondly, again motivated by ( 1271) , we would like to determine the dependence of the 
coherence lengths £ on the chemical potential fi. This is shown in Fig. [TJJ It is seen 
that for values of the chemical potential close to the critical value (i.e., T close to T c ), 
these mass scales are linear functions of the chemical potential (in this manuscript, 
we use the term mass scale interchangeably with inverse length scale). This is true for 
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Figure 12: The coherence lengths as a function of the chemical potential (0\ on the 
left and O2 on the right) 



both the 0\ and 2 solutions - but £ and £ g have slightly different slopes as functions 
of the chemical potential. 

This is in sharp contrast to what one might have expected from either GP-equation 
or even the bulk equations. The numerical solutions display different length scales 
for the condensate and charge density (for fi > fi c ). We have confirmed that this is 
not a numerical artefact in the following ways. As a first check, we plot the behavior 
of the two length scales as we increase lattice size in Fig. [T3J We see that the 




Figure 13: The best fit valued for £ (bottom) and £ g (top) vs. lattice size for an C 2 
condensate at — = 1.03. 

difference between the £'s saturates for large lattice sizes, indicating that errors due 
to discretization will not swamp the difference between the two scales. 

We also plot the relative difference between £ and £ g as a function of the chemical 
potential in Fig. [TU Following the earlier discussion on error estimate for the con- 
densate and charge density profiles, we have determined error bars on these numbers 
(which have also been indicated in Fig. UM . Since the error bars do not overlap with 
zero (except close to T c ), it is clear that the difference in these length scales is not a 
numerical artefact. A sharp difference between the two kinds of condensates is seen 
in the observation that the length scale difference changes sign. 
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The existence of two distinct length scales leads one to expect an interesting 
layering effect as one nears the soliton's core. Two length scales would indicate 
that the relative fraction of charge density in the condensate versus non-condensate 
degrees of freedom varies as a function of x. For 0\ we find that the condensate's 
length scale is larger than that of the total charge density. As one approaches the 
soliton from infinity, the 0\ condensate starts to drops off before the total charge 
density, therefore the non-condensate contribution to the total charge density must 
be must be relatively over dense when compared to its asymptotic value. For O2, 
the relative order of the length scales is reversed, implying that the density in the 
non-condensate fraction must be relatively under dense near the dark soliton. 

4.3 Depletion fraction 

Another quantity of interest in these objects is the amount of density depletion at 
the core of the soliton. This is illustrated in the left panel of Fig. [TS] where we 
plot the fractional density as a function of the distance from the core of the soliton 
(for -p- = 1.9). One sees that there is a striking difference in the amount of density 
depletion at the core between 0\ and 2 type of condensates. Fig. [15j In the 
right panel, we plot the percentage density depletion at the core as a function of 
the chemical potential for the two kinds of condensates. Firstly, the fact that this 
depletion fraction is not 100% is an indicator that we are quite far away from the 
solution ff26]) of the GP equation. 

In [20], it was shown that the density depletion fraction, at zero temperature, 
was directly related to whether the system was BEC-like (large depletion) or BCS- 
like (small depletion) (for non-relativistic systems). It is interesting to note that the 
density depletion for holographic solitons in Fig. [15] is strongly dependent on the type 
of operator which condenses. The amount of density depletion for the 0\ condensate 
is likely to be quite large near zero temperature. By contrast, O2 type condensate 
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0.5 0.6 0.7 0.8 0.9 1.0 T c 



Figure 15: Left: density depletion: 0\ is blue and dotted. Right: density depletion 
as a function of the temperature: 0\ on the top and 2 on the bottom 

seems to be saturating near 40% density depletion as we tend to zero temperature. 

If one defines the type of superfluid (BEC or BCS) by the amount of density 
depletion at the core [20] , then Fig. [15] suggests that the 0\ system is likely to be 
BEC type. However, as we shall see later, it is to be emphasized that the comparison 
with the BCS/BEC types of superfluids is perhaps an analogy only 

4.4 Comparison of solitons with quasiparticles 

The dark soliton we have constructed is a "macroscopic" object in this system with 
a characteristic length scale, namely the coherence length. On the other hand, we 
also have quasiparticles which are massive. These are "microscopic" excitations, and 
one could try to compare the inverse mass of the quasi particle with the coherence 
length. In a sense, since linear response is controlled by the quasiparticle mass, it 
will be interesting to wonder if the coherence length may be accounted for using some 
heuristic based on linear response theory. 

The quasiparticle mass can be determined by studying the two point functions 
for fluctuations around the homogeneous condensate, and examining the fall-off as 
a function of the spatial co-ordinate x (we define the quasi-particle mass this way). 
The boundary conditions satisfied by the fluctuations in this case are "reflective" 
regularity conditions at the horizon. This is different from the sound mode studies in 
that the latter uses infalling boundary at the horizon as is appropriate for a black hole 
quasi- normal mode (which represents a relaxation mode for the superfluid system). 

Using the AdS / CFT dictionary, the boundary quasiparticle masses arise from the 
poles in the bulk to boundary propagator in Euclidean AdS space [21]. Because we 
want the lightest masses (largest length scales), it is only necessary to solve for the 
lowest pole which has vanishing Matsubara frequency. Finding the lowest pole of the 
bulk to boundary propagator reduces to finding the values of k 2 for which one has a 
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static solution to the linearized equation of motion satisfying the correct boundary 
conditions (see below). 

In detail, we can linearize ([7j) to find 

A 2 2 AH 

fSR" + f'5R' - z5R - k 2 5R + 6R{--) + ~J~ 5A = ( 28 ) 
f6A" - k 2 5A - R 2 5A - 2AR5A = (29) 

where k represents the (dimensionless) momentum of the particle in the x-direction 
and we have dropped the tilde on the fluctuation SR(k, z) (R, A are solutions rep- 
resenting the background). The presence of the interation term between 5R and S A 
suggests that there is mixing between these two kinds of quasiparticles - very similar 
to the mixing between axial and vector-mesons in models of AdS/QCD [22]. 

The boundary conditions at the horizon are again obtained by requiring regularity 
of the solutions 

- 35R'(k, z = 1) = (1 + k 2 ) 5R(k } z = 1) 5A(k,z = 1) = (30) 

One may understand this as a requirement that the excitation energy of the quasi- 
particle over the ground state energy be finite. 

In order to determine the boundary conditions to be imposed at the AdS boundary, 
we again determine the behavior of the modes close to the boundary and require that 
the corresponding "non-normalizable" modes are zero. This gives 

5R( Z = 0) = 8A(z = 0) = (31) 

We now solve the equations (1291 subject to the above boundary condition as an 
eigenvalue problem for M 2 = —k 2 . For this value of k 2 there will be a vanishing 
eigenvalue for the linearized equations of motion, and hence a pole in the bulk to 
boundary propagator. The resultant graph of quasiparticle mass vs. temperature is 
shown in Fig. [TH] along with the coherence lengths. It is interesting to note that for (9 2 
the quasiparticle length scale tracks the soliton's condensate length scale for — < 1.2. 
This is surprising in that one might have expected the lightest quasiparticle to follow 
the largest length scale. This might be a sign that one should also include other fields 
in the quasiparticle analysis (in the sense that one of the other fields has poles which 
track the larger length scale). To date, our numerical accuracy is insufficient to see 
if this is also true for 0\ . 

5 Discussion 

In this work, we have studied properties of one type of extended configuration allowed 
in the holographic superfluids described in [6]. These superfluids support dark soliton 
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Figure 16: Quasiparticle mass (magenta, dot-dashed) as a function of fi with £ (blue) 
and £ q (red, dashed) 

: 0\ on the left and O2 on the right. 

solutions which are characterized by a local depletion in the charge density. This paper 
explored various features of these solitons as functions of the chemical potential and 
condensate type using numerical simulations. 

Although the system we actually solved is a Maxwell-Higgs system in the presence 
of a black hole, the magic of the holographic AdS/CFT correspondence results in a 
dark soliton with parallels to the soliton of the GP equation. In fact, we found that 
the length scale associated with the soliton scales with the chemical potential as for 
soliton of the GP equation. 

However, we found that the variation in the order parameter occurs with a different 
length scale than for the charge density. This is a feature for superffuids of both 
(Oi) and (O2) types. The presence of these two length scales is a surprising result 
which is not manifest in the original gravitational system nor was it predicted by 
hydrodynamical studies. 

A feature of these scales is that for the (0\) superfluid, the length scale associated 
with the order parameter is larger than that associated with the total charge density. 
As discussed in Section \A.2\ this indicates an interesting spatial dependence of the 
non-condensate fraction of the charge density. If we use the size of order parameter 
as an indicator of the fraction of charge density which resides in the condensate itself, 
we can conclude that the non-condensate fraction of the total charge density is over 
dense near the dark soliton. For the (O2) superfluid, the relative size of the two length 
scales is reversed and we find that the non-condensate fraction must be under dense 
near the soliton core. 

Spatial ordering in the components of the charge density will almost certainly 
affect transport phenomena near the dark soliton. The presence of different spatial 
orderings in the non-condensate charge density might allow one to test how the quasi- 
particles see the condensate differently than the non-condensate matter, although it 



21 



should be pointed out that the relative difference in the two length scales is small for 
the chemical potentials discussed in this paper. This would be an interesting direction 
to pursue in future work. 

A second feature that we found is that, at finite temperature, the density deple- 
tion in the core of the soliton is typically quite far from 100%. In this sense too, 
the solitons of the GP equations are a bad guide. The actual value of the density 
depletion fraction strongly depends upon the type of holographic superfluid studied. 
For an (0±) condensate we find that the depletion fraction grows to near 100% for 
low temperatures, whereas for an ((9 2 ) condensate the depletion fraction seems to 
asymptote towards 40% for low temperatures. 

These same gross features are known to exist in non-relativistic superfluids. In 
the context of the BEC-BCS crossover, it is known that BEC dark solitons have near 
100% depletion fractions, while BCS superfluids have depletion fractions less than 
60% [20]. In that case, the difference is associated with the characteristic size of the 
field which condenses. In the BEC regime it is a pointlike boson condensing while in 
the BCS regime the condensate is comprised of large Cooper pairs. 

It should be noted that in holographic superfluids, the precise microscopic de- 
scription of the system is not known. In general, it may be expected that we can 
construct both fermionic and bosonic quasiparticles. One could ask if an analogous 
classification of the size of the quasiparticles (into tightly bound BEC atoms and 
relatively large sized Cooper pairs) is possible in relativistic holographic superfluids. 
Therefore, we could hope to realize features of both types of both BEC and BCS type 
superfluids by tuning parameters of the gravitational description. 

If we associate the depletion fraction of dark soliton with the scaling dimension of 
the condensing operator, it is natural to imagine that one may change the depletion 
fraction as one tunes the scalar mass in the gravitational action. In fact, in the works 
of [16,23], the authors found a rich structure of ground states as the mass of the bulk 
scalar varied away from L 2 m 2 = —2. It would be quite interesting to explore how the 
properties of the quasiparticles and solitons vary with m 2 and whether one can find 
a crossover analogous to what is found for non-relativisitc superfluids. 

Of course, one of the important question to be answered before taking the analogy 
to the BCS system seriously is whether there is a Fermi surface, or some aspect of 
a Fermi surface in the dual system. There have been several explorations in this 
regard and in fact, it has even been studied if one can realize a (non-)Fermi liquid 
holographically [24-27,29,30]. 

One potential shortcoming of our work is the fact that we do not have a smooth 
limit to zero temperature. This is an artefact of the probe limit [31]. Therefore, a 
numerically challenging problem which is potentially of much interest is to construct 
the fully backreacted black hole with a dark soliton. Such a solution to the gravity 
system would be very interesting in the sense that it would give this black hole rather 
novel hair. For a recent discussion along this line see [32] (of course, in string theory 



22 



contexts, such spatially inhomogeneous hair have a complex history). 

We have observed (I2.2p that the presence of the condensate perturbs the gauge 
field background by little. Therefore, it might be interesting to explore a different per- 
tubation expansion, by considering scalar fields in a charged black hole background. 
This would amount to including the backreaction of the noncondensate part, and thus 
correspond to an expansion in the superfluid density. 

We also note that in recent works [33-35], several other extended solutions of 
the gravitational system were studied. In these works, non-normalizable components 
of some bulk fields were turned on and the solitons were then interpreted as being 
vortices in a superconductor. In this sense, while these authors study the same set 
of equations, their results are interpreted quite differently. In a subsequent work, 
we will present numerical results about vortices in a supefluid. The main difference 
with the previous work is the absence of any non-normalizable component to the bulk 
magnetic field. 

Finally, one would like to explore such phenomena in the context of non-relativistic 
gauge gravity duals. In this way one might hope to model the experimentally observed 
dark solitons much more closely. 
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